
# Plot clone size distrib over time
par(mfrow=c(1,1))
for(i in 1:1){
        sizedistr <- read.csv(paste("sweep",i,"/sweep",i,".unique.clones.size.distrib.log",sep=""),header=FALSE)
	simend <- sizedistr[sizedistr[,1]==7300,2]

	# Exclude size one
	simend <- simend[simend>1]

	xlimit <- 100

	hist(simend, freq=FALSE, main="Clone Size Distribution",breaks=600,xlim=c(0,xlimit),xlab="Clone Size",col="gray",border="black")
	distmean <- mean(log(simend))
	distsd <- sd(log(simend))
	curve(dlnorm(x,mean=distmean,sd=distsd),from=0,to=xlimit,add=TRUE,col="red",lwd=2)
	distmean <- mean(simend)
	distsd <- sd(simend)
	curve(dlnorm(x,mean=distmean,sd=distsd),from=0,to=xlimit,add=TRUE,col="blue",lwd=2)
	rate <- mean(1/simend)
	curve(dexp(x,rate=rate),from=0,to=xlimit,add=TRUE,col="darkgreen",lwd=2)

}


for(i in 1:1){
        distribs <- read.csv(paste("sweep",i,"/sweep",i,".unique.clones.size.distrib.log",sep=""))
	# A clone is a unique ms pattern
        clones <- read.csv(paste("sweep",i,"/sweep",i,".grid",sep=""))
	# pick a random color for each clone
	clonecolors <- NULL
	clonecount <- max(clones[,2])
	red <- sample(c(0:255),clonecount,TRUE)
	green <- sample(c(0:255),clonecount,TRUE)
	blue <- sample(c(0:255),clonecount,TRUE)
	clonecolors <- cbind(c(1:clonecount),red,green,blue)
	# plot the matrix
	mat <- matrix(crypts[,1],nrow=256,ncol=256,byrow=TRUE)
	matcol <- mat
	for(x in 1:256){
	for(y in 1:256){
	if(x==256&y==256){
		matcol[x,y]<-rgb(255,255,255,max=255)
	} else {
matcol[x,y]<-rgb(clonecolors[clones[(x-1)*256+y,2],2],clonecolors[clones[(x-1)*256+y,2],3],clonecolors[clones[(x-1)*256+y,2],4],max=255)
	}
	}
	}
par(mfrow=c(1,1))
par(pty="s")
color2D.matplot(mat,cellcolors=matcol,xlab=NA,ylab=NA,do.hex=TRUE,border=NA)
}



for(i in 2:2){
        sizedistr <- read.csv(paste("sweep",i,"/sweep",i,".unique.clones.size.distrib.log",sep=""))
        boxplot(split(sizedistr[,2], sizedistr[,1]),ylim=c(0,2000))
        clonenums <- as.numeric(sapply(split(sizedistr[,2], sizedistr[,1]),length))
}

library(plotrix)

# Plots 5x5 hexagonal grid
mat <- matrix(c(1:25), nrow=5, byrow=TRUE)
back <- matrix(rep("gray",25), nrow=5, byrow=TRUE)
# Make square plotting region for even hexagons
par(pty="s")
color2D.matplot(mat,cellcolors=back,xlab=NA,ylab=NA,do.hex=TRUE)

# Plot 
par(mfrow=c(2,2))
mat <- matrix(c(1:25), nrow=5, byrow=TRUE)
back <- matrix(rep("gray",25), nrow=5, byrow=TRUE)
deadcells <- back
deadcells[3,2] <- "white"
deadcells[4,2] <- "white"
deadcells[3,3] <- "red"
color2D.matplot(mat,cellcolors=deadcells,xlab=NA,ylab=NA,do.hex=TRUE,axes=FALSE)
deadcells[4,2] <- "red"
color2D.matplot(mat,cellcolors=deadcells,xlab=NA,ylab=NA,do.hex=TRUE,axes=FALSE)
deadcells <- back
deadcells[3,3] <- "red"
color2D.matplot(mat,cellcolors=deadcells,xlab=NA,ylab=NA,do.hex=TRUE,axes=FALSE)
deadcells[2,3] <- "red"
color2D.matplot(mat,cellcolors=deadcells,xlab=NA,ylab=NA,do.hex=TRUE,axes=FALSE)


# Plots 64x64 hexagonal grid
mat <- matrix(c(1:(32*32)), nrow=32, byrow=TRUE)
back <- matrix(rep("white",32*32), nrow=32, byrow=TRUE)
# Make square plotting region for even hexagons
par(pty="s")
par(mfrow=c(1,1))
color2D.matplot(mat,cellcolors=back,xlab=NA,ylab=NA,do.hex=TRUE)


par(mfrow=c(2,1))
hist(rexp(1000,1/4.0),xlim=c(0,25),ylim=c(0,0.25),breaks=15,main="Time to Cell Division",freq=FALSE)
hist(rexp(1000,1/6.0),xlim=c(0,25),ylim=c(0,0.25),breaks=20,main="Time to Cell Death",freq=FALSE)

par(mfrow=c(1,1))
curve(dexp(x,1/4.0),from=0,to=20,col="red")
curve(dexp(x,1/6.0),from=0,to=20,add=TRUE,col="blue")



